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GAP OPENING BY EXTREMELY LOW-MASS PLANETS IN A VISCOUS DISK 
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ABSTRACT 

By numerically integrating the compressible Navier-Stokes equations in two dimensions, we calculate 
the criterion for gap formation by a very low mass (q ~ 1CP 4 ) protoplanet on a fixed orbit in a thin 
viscous disk. In contrast with some previously proposed gap-opening criteria, we find that a planet 
can open a gap even if the Hill radius is smaller than the disk scale height. Moreover, in the low- 
viscosity limit, we find no minimum mass necessary to open a gap for a planet held on a fixed orbit. 
In particular, a Neptune-mass planet will open a gap in a minimum mass solar nebula with suitably 
low viscosity (a < 1CP 4 ). We find that the mass threshold scales as the square root of viscosity in 
the low mass regime. This is because the gap width for critical planet masses in this regime is a fixed 
multiple of the scale height, not of the Hill radius of the planet. 

Subject headings: hydrodynamics - methods: numerical - planet-disk interactions - planets and sat- 
tellites: formation - protoplanetary disks 



1. INTRODUCTION 

The question of whether a gap will form in a planet- 
disk system i s of critical importan ce in understanding 
its evolution (jKley Sz Nelson ll2012| ) . The amount of gas 
in the planet's immediate vicinity will st rongly deter- 
mine the migration rate through the disk dWard Ill997f) 
and potentially even halt migration dWard fc Hourigan 
[19891 ILi et al. I [20091 lYu et al. I Mm . Gap opening 
can also dramatically affe ct accretion onto the secondary 
(|Lin fe Papaloizo"u~lll993ft . As such, planetary evolution 
models are sensitive to the gap opening criterion, and it 
is therefore of great interest to know this criterion pre- 
cisely. 

Se veral conditions have bee n proposed for gap open- 
ing. iLin fc Papaloizoul (|1993ft proposed a stability limit, 
which is essentially a limit on the shear produced by 
steep pressure gradients in the disk. The Rayliegh Sta- 
bility criterion reduces to a condition that the planet's 
Hill radius be larger than the disk scale height: 

R H = K9/3) 1 / 3 > h, (1) 

where r is the orbital radius, q = M v /M* is the mass 
ratio, and h is the disk scale height. We will refer to 
this limit as the "Strong Shock" limit. This is because it 
is related to the requirement that a strong shock forms 
within a scale height of the planet's orbit. In terms of a 
limit on the planet mass, ([T]) reduces to 



M p > 3M Sh , 
M Sh = hc 2 /G = M*/M £ 



(2) 
(3) 



Here, M p is the mass of the secondary, c is the local 
sound speed, M* is the primary mass, and Ai = r/h 
is the disk Mach number. ILin fe Papaloizou (1993j) also 
demonstrated a viscous limit, which requires that the gap 
be opened faster than viscosity can refill it: 



q > 



40v 



(4) 



where q is the mass ratio, v is the viscosity, r is the dis- 
tance between the primary and the secondary, and fi is 
the orbital frequen cy. In terms of the Shakur a-Sunyaev 
viscosity v — ach (Sha kura fe Sunvaev l Fl973). this vis- 
cous limit is 

M p > M a , (5) 



AOaMMsh- 



(6) 



Another important gap opening criterion is t he inertia! 
limit , which appropria te for migrating planets (|Li et al. I 
2009: lYu et all 12010). In this case, the gap opening 
rate must be faster than t he secondary's migration rate. 
iWard fc~H ourigan (1989) calculated this limit: 



M p > Mi, 
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(7) 
(8) 



Here, S is the surface density. In typical protoplanetary 
disks, the disk mass Er 2 is much smaller than the pri- 
mary mass, meaning that planets satisfying the strong 
shock criterion are essenti ally guaranteed to s atisfy ([5]). 
Synthesizing these results, iCrida et al. I (|2006[ ) built an 
effective ID description of the disk profile and used it 
to construct a gap opening criterion which combines the 
viscous limit with the strong shock limit: 
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Note that already the "strong shock limit" has been re- 
laxed to M p > l.SMsh- However, this picture cannot 
be acc urate. It has been demonstrated, both the oret- 
ically (|Goodman fc Rafikov I [20011 iRafikov 1 12002a!) and 
in dire c t numerical calculations {Dong, Rafiko v fc Stone"! 
I2011H iDuffell fc MacFadven 1 12012ft that planets with 
masses well below Msh ca n still op e n gaps . 

The point, as noted by iRafikovl (|2002bft . is that gap 
opening does require a dissipation mechanism such as 
a shock, but not necessarily a strong shock. In fact, 
any perturbing mass on a fixed orbit in an inviscid disk 
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will eventually open a gap. However, in real situations 
with nonzero viscosity, "eventually" may entail unrealis- 
tic timescales, i.e. longer than the viscous timescale or 
the migration timescale. 

If the strong shock limit is not a necessary condition for 
gap opening, then migration effects become important, 
and the inertial limit © would become relevant. Ad- 
ditionally, three-dimensional effects could no longer be 
neglected, because the planet's radius of influence would 
not extend across the entire thickness of the disk. More- 
over, the viscous limit (|6|) was derived assuming that the 
strong shock criterion was satisfied, so it might take a 
different form when applied to low mass secondaries. 

To understand why this might be the case, we can es- 
timate the rate of angular momentum transfer due to 
the planet a nd due to viscosity, as ha s been previously 
estimated bv iLin fc Papaloizou I (|1993| ): 

*--*(*)'(£)' <-» 

H a ~ vllr 2 ~ achZr 2 , (II) 

where A is the gap width, indicating that the tidal 
torque depends on how far the edge of the gap is from 
the secondary (This crude estimate can be improved 
upon by direct calculations of planet torque, see for ex- 
ample lETAngek) fe Lubow I (pOlO ). iRafikov fc Petrovich I 
(2012)). The gap opening criterion can be estimated by 

equating these torques, H p ~ H a : 

If the gap width is of order the disk scale height A ~ h, 

M gap ~ MshV^M (13) 

However, if the gap width is governed by the secondary's 
Hill radius, A ~ Rh ~ rq 1/3 , 

M gap ~ M sh {aM), (14) 

which is in agreement with the form f ound for large mass 
planets bv ILin fe Papaloizou I (|l993), our equation ©. 
This indicates that the scaling with a is dependent upon 
whether the width of the gap is governed by the scale 
height of the disk, or the Hill radius of the planet. One 
might suspect that it is governed by whichever of the 
two is greater, meaning that if the secondary is small 
enough that its Hill radius is smaller than a scale height, 
the width of the gap is independent of the mass, and the 
scaling should correspond to (fl"3|) . This turns out to be 
correct, as we shall confirm in this work. 

Using a direct numerical approach, we find the vis- 
cous limit for a low-mass secondary (one which does not 
satisfy the strong shock criterion). We keep the planet 
at a fixed orbital radius and vary its mass and the disk 
viscosity, to empirically determine the gap opening crite- 
rion as a function of the dimensionless Shakura-Sunyaev 
viscosity parameter, a. Migration effects are intention- 
ally neglected, so that we may isolate the dependence on 
viscosity alone. 

The numerical inte gration is carried out using a movin g 
computational mesh (|Duffell fc MacFadven ||20111|2012D . 
This reduces numerical viscosity and allows us to take 



long time-steps, thus making it possible to run the calcu- 
lations for a large number of orbits, and hence to explore 
the low- viscosity regions of parameter space. 

2. NUMERICAL INGREDIENTS 

The equations being solved are the compressible 
Navier-Stokes equations. For an effective 2-D descrip- 
tion of the gas, these equations are vertically integrated 
(e.g. £ = J pdz). In other words, vertical structure is ne- 
glected, and vertically propagating modes are considered 
subdominant: 

dt^ + d^v,) = (15) 

dtptoj) + dipViVj) + d 3 P = FJ lsc + Ff av (16) 

d t E + di((E + P) Vi ) = (Fp sc + F° rav ) Vj (17) 

Where £ is the surface density, v is the velocity, P is the 
pressure, and E is the energy density, E = l/2pv 2 + e int . 
The equation of state is taken to be adiabatic, 

P = ei„t (7-1) (18) 

but we choose 7 = 1.001 to make the equation of state 
nearly isothermal. 
The source terms include the viscous force, 

Fp sc = diip^vidiVj + djVi) + CSijdkVk)), (19) 

where v is the kinematic viscosity, and the gravita- 
tional force, F grav , which is the force produced by both 
point masses. These masses are each moved in a Keple- 
rian circular orbit about a common center of mass. In 
our implementation, the viscous force is re-expressed as 
a viscous flux and moved to the left hand side of the 
equation. 

2.1. DISCO Code 

DISCO is a finite-volume, moving-mesh hydrodynam- 
ics code which is specifically tailored to the study of 
gaseous disks. DISCO utilizes a cylindrical grid which 
moves and shears azimuthally with the orbital velocity 
of the fluid. It is based on the TESS code, which is a gen- 
eral moving-mesh method for solvin g hyperbolic systems 
(|Duffell fc MacFadven lf2?TTlL l20"H . 

The computational domain runs from r — 0.4 to r — 
1.6 with 512 radial zones, which is about 25 zones per 
scale height in the vicinity of the planet's orbit. The 
number of azimuthal zones varies with radius, in order 
to keep a fixed aspect ratio of 1:1. This amounts to 
« 3200 zones. 

The calculation of the viscous terms is detailed in the 
appendix. Visco sity on a moving me sh has already been 
implemented by iMunoz et al I ((2012) for a Voronoi tes- 
sellation, but because our time integration utilizes the 
method of lines, our implementation is much simpler. 

Kinematic viscosity in the disk is assumed to be a con- 
stant, but we report our results in terms of the dimen- 
sionless parameter a = vc p h p (the subscript p indicates 
we are evaluating these quantities at the planet's orbital 
radius). In other words, a is not spatially constant, but 
we report our results in terms of its value in the vicinity 
of the secondary. 
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2.2. Background Flow 

The initial cond itions assume th e Minimum Mass Solar 
Nebula (MMSN) (jHavashi 111981 : 



£(r) = So(r p /r) 3 / 2 
c(r) = c (r p /ry/ 4 



(20) 
(21) 



with orbital frequency set to balance pressure and grav- 
itational forces: 

n 2 (r)=nt(r)(l-2^r~/r~ p /M 2 y (22) 



fio(r) = VGA/*A 3 



(23) 



Because we keep the planet at a fixed orbit and the gas 
exerts no gravitational force, So is arbitrary, cq is fixed 
by the mach number: 



c = r p tt /M. 



(24) 



For the Mach number, we choose Ai = 20, which cor- 
responds to an orbital radius r p = 2.4 AU. The Mach 
number is not fixed as a function of orbital radius. 

2.3. Perturbing Potential 

The gravitational well of the planet is given by the 
simple Newtonian formula: 



$(s) = GM p /s 



(25) 



where s is the distance from the planet to the fluid el- 
ement. However, we use a smoothed potential, for two 
reasons. First, the formula (|25l) is divergent, and its sin- 
gularity is on the grid, which would make time integra- 
tion impossible. Secondly, in this effective 2-D treat- 
ment, (|2~51 would overestimate the secondary's influence 
within a scale height of the planet, because our two di- 
mensional fluid elements actually feel a vertical average 
of this gravitational f orce. To account for this, as has 
been done by others dTanaka. Takeuchi fe Ward 1 120021 : 
iMasset I l2002t iMuller. Kiev fe Meru I l2012f >. we use an 
approximate "vertically averaged" potential of the form 



$(a) = GM p /y/s 2 + S 2 



(26) 



where we choose S = .6h. The inner and outer radial 
boundaries are ha ndled via the same damping p rocedure 
that was used in IDuffell fe MacFadvenl (|2012t ). As in 
that work, we do not attempt to model accretion onto 
the planet (The effects of accretio n on the gap forma- 
tion p rocess have been discussed bv lD'Angelo fe Lubow I 

3. RESULTS 

3.1. Saturated Gap Depth 

Each system evolves for thousands of orbits before fi- 
nally saturating, the secondary having hollowed out a 
gap with some minimum density. We plot minimum 
density on the grid as a function of time for the AMsh 
secondary in Figure [TJ for various viscosities. For this 
massive secondary, saturation is established after a few 
thousand orbits, and the saturated minimum is propor- 
tional to the viscosity. We plot the final E in all test 




500 1000 1500 2000 2500 

Time (Orbits) 

Fig. 1. — After many orbits, the planet opens a gap, which 
eventually saturates to its final steady-state depth, at which the 
rate of viscous filling is equal to the rate of evacuation. Plotted here 
is the (time-averaged) minimum density in the gap, as a function 
of time, for various disk viscosities. The final saturated minimum 
is proportional to the viscosity. The planet mass in this example 



is M„ 



AMg h (About half of Jupiter's mass). 



cases in Figure [5] These can be mapped to simple scal- 
ing relations, as a function of M p /Msh and a. The min- 
imum density is linear in a and inversely quadratic in 
the secondary mass, as it is a nonlinear effect of the per- 
turbation. We construct an empirical scaling relation for 



(27) 



We note that the gap depth obeys power-law scaling 
with respect to viscosity and planet mass. This result is a 
bit unexpected from th e standpoint of ID torque-balance 
models of gap profi les (IVarniere. Quillen fe Frank~ll2004t 
iCrida et al. I l2006h which typically predict exponential 
dependence of gap depth on system parameters. 

Equation ((27)) should be slightly modified to take into 
account secondary masses below the threshold, for which 
this formula would erroneously predict T, sat > So- Re- 
placing T, sat with 



Ssat — SsatSo/ (Esat + Ed) 



(28) 



fits the low-mass, high- viscosity data reasonably well 
(Figure^]). However, the behavior for below-critical per- 
turbers is less important; equation (|2"7) gives us the ma- 
chinery to address the question of gap opening. 

3.2. Gap opening criterion 

We choo se our diagno s tic cr iterion for gap opening to 
agree with lCrida et al. I (|2006| ): 



Ssat < O.lEo, 

Using (|2~T|) we arrive at the criterion 



(29) 



M gap = M S h(17VaM). 



(30) 

We plot this criterion in FigureEJ alongside data points 
which were interpolated from Figure [5J We find ex- 
cellent agreement. Alongside this, we plot th e "Strong 
Shock" prediction of lLin fc Papaloizou l (119931) . and the 
combined prediction of ICrida et al" J d2006() (Equation 
[9]). We note what appears to be disagreement between 
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Fig. 2. — The saturated gap depth obeys simple scaling formulas, 
as a function of disk viscosity a and the mass of the secondary. We 
plot minimum gap density for all test cases. The dashed curves 
represent the sca ling relation 127t . Solid curves incorporate the 
modification (1 28 6 - 



Gap Opening Threshold at 2.7 AU in the MMSN 
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Fig. 3. — Gap opening criteria, dependent on viscosity and planet 
mass. We plot several curves representing potential gap opening 
criteria. Planet-Disk systems below a given curve will not open 
a gap according to the criterion. We plot equations J2I and (0, 
corresponding to predicted gap-open ing thresholds, alongside our 
empirically derived scaling relation 1 1301) . To highlight the accu- 
racy of our scaling relation, we plot two distinct measured criteria 
for gap-opening. Open squares represent the viscosity at which 
the minimum surface density dips below £ < O.lErj, an d the filled 
squares represent the viscosity at which the total torque on the 
planet changes sign, a noisier measurement which is nevertheless 
roughly consistent with the minimum surface density measurement. 



these previous predictions and our current empirical re- 
sults, at least for low viscosities and small planets. Eq. 
^ asymptotes to a finite value in the low-viscosity 
limit, corresponding to the "Strong Shock" limit M p > 
l.lMsh- This discrepancy could be numerical in ori- 
gin; an under-resolved calculation with a small but non- 
negligible numerical viscosity (a num ~3x f would 
spuriously produce this behavior. This emphasizes the 
importance of low numerical viscosity when studying this 
low-mass planet regime. 

The scaling given in (|50|) c oincides with older esti- 
mates for high mass pla nets ()Lin fc Papaloizoul 119861 : 
IWard fc Hourigan1ll989D . The idea that this same scal- 
ing would apply for low-mass planets was suggested by 
iRafikov 1 (|2002b| ) . It should be noted that in this work we 
have restricted ourselves to the case M — 20; the scaling 
with Mach number should be checked in a future work. 
This formula suggests the scaling (fl"3|) . which implies a 
fixed gap width A ~ h. This can be investigated by ex- 



Gap Profile for a M p = .5 M Sh Secondary 
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Fig. 4. — Azimuthally averaged surface density at late times 
for a planet of roughly Neptune's mass (M p = .5Msh)- Curves 
are plotted with various viscosities. Thou gh the planet mass is 
well below the limit of Crida et al. (2006), a gap opens at low 
viscosity (a < 10~ 4 ). 

amining the gap profile directly. 

3.3. Gap profile 

We plot the azimuthally averaged density profile for 
one planet, M p — .5Msh, with various viscosities, in Fig- 
ure [U Though it is too small to satisfy the strong shock 
condition, clearly this secondary opens a gap when the 
viscosity is low enough. 

The scaling relation (f30|) . when compared with (|13|) 
and pT|) , seems to imply for low- mass planets that the 
gap width should be given by the scale height , and not 
the Hill radius of the secondary. In order to study this, 
we plot the gap profile in the low- viscosity limit for our 
various planet masses (Figure [S]). Though the largest 
mass is 16 times the smallest mass, there is not a signif- 
icant difference in gap width between the different pro- 
files. However, as there is some small variation in gap 
width, we plot its dependence on planet mass in Figure 
El Here, the edge of the gap is defined as the radial posi- 
tion at which £ < E /3. The width is plotted in Figure 
[5] for various masses and viscosities, and it is clear there 
is a weak dependence on planet mass. In the right panel 
of this figure, each curve is re-scaled so that the planet 
mass is expressed as a fraction of M gap (Equation [30]) . 
Interestingly, the curves appear to converge at the criti- 
cal mass M p = M gap , at a fixed width of A ~ 6h. This 
suggests that the shape of the gap profile at the critical 
mass may be r oughly in d epende nt of other parameters 
of the system. IRafikov 1 ()2002al ) suggests that the gap 
edge should coincide with wherever the density wave dis- 
sipates (where it shocks, in the inviscid case). 

3.4. Stability of Deep Gaps 

Looking again at Figure [TJ we note that the inviscid 
calculation is plotted alongside the viscous results, ap- 
parently having a slightly lower effective viscosity than 
the a = 3 x 10~ J run. Naively, we might be tempted 
to interpret this effective viscosity as numerical, but we 
shall see in the next section that our numerical viscosity 
is ctnum ~ 10~ 6 , which is much too small for this inter- 
pretation. 

A hint as to the explanation can be seen in Figure 
Q] at around 1200 orbits. The inviscid calculation has 
some intermittent, erratic behavior which settles down 
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Fig. 5. — Azimuthally averaged surface density at late times for 
various planet masses (a = 3 X 10 ). For low mass planets, the 
gap width is weakly dependent on planet mass. 
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Fig. 6. — The gap width is given for various planet masses and 
viscosities. The edge of the gap is defined as the point where £ = 
So/3. In the right panel, we re-scale the pl anet mass as a fraction 
of the critical gap-opening mass (Equation I30D ■ We find that the 
curves for various viscosities converge at around M p ~ M gap , with 
a fixed width of A ~ 6h. 

on roughly a 10 orbit timescale. This short-timescale 
jostling occurs about every thousand orbits or so, and it 
is due to vortices forming and breaking up after several 
orbits. These vortices seem to provide an effective viscos- 
ity which appears to be the cause of the saturation floor; 
When we double the resolution this saturation floor does 
not change. Similar vortices have been witn essed and 
analyzed in previous planet-disk calculations (|Li et al. I 
l2009tlYu et al. ||2010D . In particular, these works demon- 
strated the impact such vortices can have on planet mi- 
gration. 
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Fig. 7. — Inviscid calculations, showing dependence of gap depth 
on resolution (M p = 4Msh)- For Nr = 32 and 64, numerical vis- 
cosity dominates. For higher resolutions, a saturated gap depth is 
shown which is not caused by explicit viscosity or numerical viscos- 
ity. It appears to be caused by a fluid instability which generates 
vortices that drift and dissipate, preventing further gap deepening. 



instability appears to re-activate every thousand orbits 
or so. 

Results from the Nr = 32, 64 and 128 runs can be used 
to construct a formula for saturated density as a func- 
tion of the resolution. This yields a formula for effective 
numerical viscosity by comparing our result with Eq. 1271 
The effective viscosity we find is 

a num = 2.5 x 1(T 3 (31) 

where Ar measures the radial size of a zone near the 
orbital radius. 

Our numerical viscosity for the Nr = 512 runs is there- 
fore approximately «5i2 ~4x 10~ 6 , roughly an order of 
magnitude below the physical viscosity we introduce, and 
also smaller than the effective viscosity produced by vor- 
tices through shear instabilities in the inviscid case. We 
should note, of course, that this measure of numerical 
viscosity is not only problem-dependent, but diagnostic- 
dependent. That is, it is expected that some regions of 
the computational domain have more numerical viscos- 
ity than others, partly because of the variable resolution 
but partly also because the numerical viscosity is proba- 
bly larger near shocks (the scaling of numerical viscosity 
would also be linear near shocks). That caveat aside, this 
diagnostic is a reasonable measure of "effective viscosity 
for this problem" . 



3.5. Resolution Study 

To determine our effective numerical viscosity for this 
problem, we performed several inviscid tests for the 
M p — iMsh secondary at lower resolution, comparing 
32, 64, 128, and 256 radial zones to our fiducial 512 ra- 
dial zone calculation. The time history of gap depth for 
these resolutions is plotted in Figure [7] We compare this 
with Figure [1] For the cases Nr = 32 and 64, the satu- 
rated gap depth is clearly a result of numerical viscosity. 
However, for higher resolutions, the solution converges 
on a gap depth. We discussed this convergence in i)3.41 
it appears to be the result of unstable vortices provid- 
ing an effective viscosity. In Figure[7]the 512 calculation 
clearly has intermittent oscillations, indicating that this 



4. SUMMARY 

We have performed 2D numerical calculations which 
directly demo nstrate that the "s t rong shock limit" ([T]) 
proposed by iLin fc Papaloizoul (|1993l ). and used by 
ICrida et al. I ()2006[ ) is not a necessary condition for gap- 
open ing. While o u r resu lts show reasonable agreement 
with ICrida et al. I (|2006l ) for large viscosity and planet 
mass, we do not agree for low viscosities and masses. 
In particular, Neptune-mass planets are capable of gap 
opening for a < 10 ~ 4 , which is not allowed according to 
these previously proposed criteria. We also find that the 
critical mass scales as the square root of viscosity when 
applied to low-mass planets. Generally, we do not find 
any minimum mass limit for the case a — > 0. 
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Because gap opening is still possible when the Hill ra- 
dius is smaller than a scale height, we expect that 3D 
calculations will be necessary to capture disk-planet dy- 
namics in this low-mass regime. Since we have found no 
"strong shock" limit for gap opening, this means that the 
inertial limit (Equation [5]) is important in low- viscosity 
disks. Up to now, we have considered the case where 
£ was negligible in the sense that the torque on the 
planet did not change its orbit; we kept both primary 
and secondary on fixed circular orbits about the center 
of mass. In order to determine the inertial limit in the 
low-mass regime, we must relax this assumption and al- 
low the planetary orbit to respond to the torque exerted 



by the disk. All of this is planned as future work. 
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APPENDIX 

VISCOSITY IMPLEMENTATION IN DISCO 

The viscous terms in DISCO are calculated as viscous fluxes, which requires calculating the divergence of the viscous 
stress tensor in cylindrical coordinates. DISCO automatically takes into account the volume element in cylindrical 
coordinates, because we use exact geometry to calculate the volumes of cells and areas of faces. Therefore, we merely 
need to re- write the viscous stresses T^- in terms of the "effective stress tensor" T e '* , defined by: 

(v ■ T)j = -L dr (^g T ;V) + -Lv^^Tjf ) + sj (Ai) 

The components of T e ' * are the viscous fluxes as they appear in the code. Note that when we express the divergence 
of the stress tensor in cylindrical coordinates, we generate a source term in the radial component of momentum, for 
the same reason that a centrifugal source term is generated when taking the divergence of the hydrodynamical stress 
tensor. This source term will have a negligible effect on disk dynamics, but is taken into account in the code for 
completeness. The components of T e ^ are 

T e 4 f = -(H(v*« r -2n) r49 , 

T e // = -(H(r 2 V r O) (A2j 

Where v r is the radial velocity and fi is the angular frequency. The first two components listed are the flux of the 
radial component of momentum, and the remaining two components are flux of angular momentum. The source term 
for the momentum in the radial direction is 



S r = -{pv)v r /r 2 . (A3) 

Following a similar path to IMufioz et al I ()2012l ) , we calculate the v isco us flux independently from the hydrodynamic 
flux (the hydro flux uses a Riemann solver). In order to evaluate (|A2[) we need to know the gradients of primitive 
variables. For this, we use the same slope-limited gradients which are used to extrapolate primitive variables from zone 
centers to faces. We extrapolate primitive variables to either side of each face, as is done for the hydro fluxes. We then 
evaluate the arithmetic mean of the variables on either side of th e fac e, and of the gradie nts, and use these averaged 
primitive variables and gradients to calculate the fluxes given in (|A2j) . The source term (|A3| is simply calculated as 
a cell-centered qu antity. 

In contrast with IMufioz et afl (|2012[ ). we do not need to extrapolate our variables in time for a half-timestep, because 
our time integration is based on the method of lines; we achieve high order by performing the full time integration via 
a series of Runge-Kutta timesteps. Each individual step can therefore be first-order in time. This simplification means 
that our viscosity implementation can be simply described by the equations (|A2I) and (| A3|) . In principle, we should 
have to take a shorter timestep if the characteristic viscous timescale for a zone is shorter than its sound-crossing time, 
but in this work we deal with viscosities small enough that we do not need to consider this timescale. 
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